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We  have  studied  the  heating  of  the  solar  chromosphere  and  corona,  and  the 
propagation  of  acoustic  waves  through  the  transition  region  between  the 
chromosphere  and  corona.  In  order  to  place  an  upper  limit  on  the  effective- 
ness of  acoustic  waves  In  heatlng,we  have  written  and  tested  a computer 
program  to  accurately  calculate  the  propagation  and  dissipation  of  verticall 
travelling  acoustic  waves  which  Includes  thermal  conduction  and  radiative 
transfer.  In  order  to  study  heating  by  other  possible  waves.  In  particular 
magneto-acoustlc-gravlty  waves 


we  have  written  and  are  now  testing  a compu 
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ter  program  which  approximately  Includes  horizontal  motions.  In 
order  to  study  the  effects  of  horizontal  Inhomogeneltles  we  have 
started  developing  a three-dimensional  fluid  dynamic  computer  program. 


We  have  also  studied  the  propagation  of  acoustic  waves  through  the 
solar  transition  region.  We  find  that  waves  with  velocity  amplitudes 
compatible  with  observations  near  the  temperature  minimum  ( <±1  km/s) 


UNCLASSIFIED 


seCUtttTV  CLASSIFICATION  OF  THIS  FASIfmiMi  Dmta  Kniaiwd; 


3. 


SOLAR  ATMOSPHERIC  DYNAMICS 

I . INTRODUCTION 

The  goal  of  this  research  Is  to  study  the  driving  of  the  solar  atmospheric 
oscillations  and  their  role  in  heating  the  chromosphere  and  corona.  Although 
many  models  of  the  solar  atmosphere  have  been  constructed  with  acoustic  waves 
as  the  energy  source,  their  ability  to  transport  sufficient  energy  has  recently 
been  questioned.  We  made  accurate  calculations  for  the  most  favorable  case 
of  vertically  propagating  acoustic  waves.  We  found  that  they  can  not  supply 
sufficient  energy  to  heat  the  corona.  Further  calculations  with  an  accurate 
treatment  of  radiative  transfer  are  necessary  to  answer  the  question  for  the 
upper  chromosphere.  A program  to  perform  such  calculations  was  prepared  and 
they  will  be  carried  out  this  coming  year.  If  acoustic  waves  can't  transport 

the  needed  energy  to  the  corona,  can  other  types  of  waves?  It  is  clear  that 

} 

coronal  heating  depends  critically  on  the  magnetic  field.  Hence  the  role  of 
magneto-acoustic-gravity  waves  should  be  studied.  Because  magnetic  waves  can 
propagate  at  a faster  speed  than  acoustic  waves  in  a strong  magnetic  field,  mag- 
netic waves  with  velocity  amplitudes  consistent  with  observations  may  transport 
the  required  energy.  We  continued  our  development  of  a "modal"  fluid  dynamic 
code  to  calculate  the  propagation  and  interactions  of  these  waves.  We  also 
began  development  of  a fully  three-dimensional  fluid  dynamic  code  that  will  in 
addi'  ^n  be  able  to  include  horizontal  inhomegenelties. 
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II . RESEARCH 

A.  Heating  of  the  Chromosphere  and  Corona 

Three  different  attacks  were  made  on  the  problem  of  how  the  chromo- 
sphere and  corona  are  heated.  First,  the  heating  by  vertically  propaga- 
ting acoustic  waves  was  investigated  including  effects  of  radiation 
transfer  and  conduction.  This  provides  an  upper  limit  to  possible 
acoustic  heating.  Second,  investigation  of  heating  by  non-vertically 
propagating  magneto-acoustic-gravity  waves  was  continued  ''^Ing  an 
approximate  treatment  of  variations  in  the  horizontal  direction.  This 
provides  a more  realistic  analysis  of  sources  of  mechanical  heating  by 
other  than  acoustic  waves  and  of  the  interaction  of  the  waves  with  the 
mean  atmospheric  structure  and  each  other.  Third,  development  of  a 
truly  3-dlmenslonal  dynamic  code  was  begun  in  collaboration  with  Richard 
Klein  (Kltt  Peak  National  Observatory)  and  Lawrence  Auer  (High  Altitude 
Observatory) . Eventually  such  a code  will  be  able  to  calculate  models 
with  horizontal  inhomogeneity.  Progress  was  made  using  all  three  approaches, 
and  work  is  still  in  progress. 

1.  Vertically  Propagating  Acoustic  Waves 

Acoustic  waves  are  generated  dynamically  by  coupling  to  the 
convective  motions  below  the  photosphere  and  also  by  overstability 
in  the  top  of  the  convective  region  (see  review  by  Stein  and  Lelbacher, 
1974,  Ando  and  Osakl,  1975).  Can  they  heat  the  chromosphere  or 
corona?  There  is  doubt  both  observatlonally  (Athay  and  White,  1978; 
Brunner  and  McWhlrter  1978)  and  theoretically  (Jordan,  1977). 

Vertically  propagating  acoustic  waves  suffer  no  refraction  and  stand 
the  best  chance  of  dissipating  their  energy  high  in  the  solar 
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atmosphere.  The  main  energy  loss  mechanisms  for  such  waves  Is 
radiative  and  conductive  damping.  A careful  treatment  of  these 
effects  Is  needed  to  assess  theoretically  the  possible  role  of 
acoustic  waves. 

For  this  purpose  a one-dlmenslonal  fluid  dynamic  code  was 
written  which  Includes  radiative  and  conductive  heating  and  cooling 
In  the  energy  equation  and  the  transfer  equation  for  radiation  In 
the  gray  - LTE  approximation.  As  In  previous  1-D  dynamic  calculations, 
the  code  Is  Lagranglan  with  the  mass  column  density  the  Independent 
variable.  The  velocity,  position  and  density  are  first  solved  for 
explicitly  from  the  equations  of  motion,  definition  of  velocity  and 
mass  conservation  (continuity)  equations  (Rlchtrayer  and  Morton,  1969). 

Then  the  energy,  Saha  and  transfer  equations  are  solved  Implicitly  and 
simultaneously  for  the  temperature,  degree  of  hydrogen  ionization,  and 
mean  radiation  Intensity  using  Newton-Raphson  iteration  (complete  lineari- 
zation) (see  Rlchtmyer  and  Morton,  1969  and  Klein,  Stein  and  Kalkofen, 

1976).  The  opacities  used  are  those  of  Robert  Kurucz  and  the  thermal 

-7  5/2 

conductivity  Is  taken  to  be  5x10  T erg/cm  K,  without  limiting  the 
flux  because  of  maximum  electron  velocity.  The  usual  treatment  of 
the  transfer  equation,  using  a finite  differenced  form  of  the  Feautrier 
equation 
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alternates  a complete  linearization  to  get  first  order  corrections 
to  variables  and  then  a formal  solution  to  obtain  the  radiation 
field  for  the  given  temperature  structure.  The  elimination  procedure 
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in  the  formal  solution  was  slightly  modified.  In  order  to  Include  very 
small  optical  depths  in  the  corona.  The  straightforward  Gaussian  elim- 
ination scheme  would  lose  all  significance  from  the  samil  term  A t J 
compared  to  the  large  term  A r It  was  therefore  necessary  to  modify 

the  method  of  formal  solution  of  the  transfer  equation  as  follows:  The 
transfer  equation  has  the  form 

\ '^1-1  ®i  '^1  ^1  “^i+l  * "°1  • ('2) 

The  usual  procedure  is  to  calculate  elimination  matri''  s 

■ 'h  * (3) 


^1+1  “ ^1  ^1-1^  » 
where  these  are  then  used  to  back  substitute  for  the 


"l  -^l+l  + "i 


At  small  optical  depth  the  radiation  field  is  nearly  constant  and 
nearly  the  identity  matrix  1.  It  is  therefore  possible  to  gain 
accuracy  by  working  with  a matrix  equal  to  the  difference  between 
and  1. 


That  is,  let 


1 - E. 


(1  - c^) 


The  relations  to  calculate  and  are  found  by  substituting  equation  7 
into  the  difference  equation  2.  The  result  is 

CBj  + - Aj  ^i-1  [b^  + A^  + Cj,  - A^  (8) 


=■1  - - [ Bi  + Ai  - Ai  ]"^  C Dl  + Aj  ] 
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The  large  terms  in  the  sum  cancelled 

analytically.  Using  this  procedure  it  is  possible  to  calculate 
to  an  optical  depth  10  where  n is  the  number  of  significant 
digets  in  the  computer  word. 

The  computer  code  was  written  and  tested  at  NCAR.  As  soon 
as  the  link  between  Michigan  State  University  and  NCAR  Is  estab- 
lished (within  the  next  few  months)  calculation  of  models  will  begin. 
These  models  will  start  with  a solar  atmosphere  in  radiative 
equilibrium.  Waves  will  be  generated  below  the  top  of  the  convec- 
tion zone  by  a piston  whose  motion  will  have  a power  spectrum 
like  that  calculated  for  waves  generated  by  the  llghthlll  mechanism 
(Stein,  1968)  plus  the  "300  second"  oscillation.  The  calculation 
will  be  allowed  to  run  until  a statistically  steady  state  atmosphere 
is  obtained.  This  atmosphere  will  have  a self-consistent  structure 
for  the  acoustic  wave  heating  - gray  LTE  radiative  losses  plus 
conduction,  included  in  the  calculation.  The  maximum  extent  to  which 
acoustic  waves  can  produce  a chromosphere  and  corona  will  then  be 
apparent.  I anticipate  that  for  the  low  chromosphere  and  possibly 
the  upper  chromosphere  the  model  atmosphere  will  accurately  repro- 
duce the  observed  solar  atmosphere.  Where  these  models  reproduce 
the  observed  solar  atmospheric  structure  they  will  be  used  to 
investigate  the  detailed  effects  of  acoustic  waves  on  line  profiles 
in  collaboration  with  Richard  Shine  (NASA  Goddard  Space  Flight  Center.) 

2.  Magneto  - Acoustic  - Gravity  (M.A.G.)  Waves 

Acoustic  waves  generated  below  the  photosphere  are  in  fact 
not  vertically  propagating,  so  they  suffer  refraction,  which  pre- 
vents them  from  reaching  the  corona  (Kopp,  1968).  Furthermore, 
they  propagate  in  the  presence  of  a magnetic  field  which  changes 

79  'M  12  Oe 


»■ 

I 

1 

t 

I 

I 

i 


■»«»-  ■ 


their  character.  Finally,  the  observed  heating  of  the  solar  at- 
mosphere clearly  depends  sensitively  on  the  magnetic  field  strength 
( Withbroe  and  Noyes,  1977).  In  order  to  study  the  propagation 
and  dissipation  of  these  M.A.G.  waves  the  development  of  an 
approximate  3-dimenslonal  "modal"  code  was  begun  in  1975  with  the 
support  of  AF  contract  F19628-75-C-0013,  The  essence  of  this 
approach  is  to  separate  all  the  variables  into  a mean  and  a 
fluctuating  part.  The  mean  part  is  uniform  on  horizor*‘al  planes 
and  the  fluctuating  part  is  expanded  in  a complet*.  set  of  horizontal 
modes  with  amplitudes  a function  of  height  and  time.  For 
simplicity  the  expansion  is  truncated  to  a single  hexagonal  hori- 
zontal mode  of  some  given  length  scale.  When  this  expansion  is 
substituted  into  the  equations  of  motion  and  transfer,  the  result 
is  a set  of  one-dimensional  partial  differential  equations  in 
z and  t for  the  mean  variables  and  fluctuating  mode  amplitudes, 
but  with  horizontal  variations  included  via  the  imposed  hexagonal 
horizontal  structure.  In  this  way  it  is  possible  to  study  non- 
longitudinal  and  non-vertlcally  propagating  waves,  but  still  include 
non-linear  interaction  of  the  wave  with  itself  in  the  vertical 
direction,  and  still  only  integrate  one-dimensional  partial  dif- 
ferential equations  using  standard  Eulerlan  techniques.  This  idea 
comes  from  Gough,  Spiegel  and  Toomre  (1975)  who  used  a similar 
method  to  study  convection. 

The  equations  were  derived  in  collaboration  with  Shyu-Hsien 
Hsieh.  They  were  coded  by  Robert  Wolff,  using  Rlchtmyers  (1969) 
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two-step  Lax-Wendroff  method  to  solve  the  difference  equations, 
with  a flux  correction  scheme  for  eliminating  numerical  Instabilities 
(Book,  Boris  and  Haln,  1975).  Testing  of  the  code  was  begun  by 
Robert  Wolff.  So  far  only  tests  of  the  one-dimensional  propagation 
terms  have  been  made,  with  mixed  results.  Piston  driven  longitu- 
dinal compression  waves  do  form  shocks  properly.  Far  from  the 
piston  the  mean  flow  develops  into  a sawtooth  pattern  with  half 
the  period  of  the  driving  piston.  This  is  because  the  piston  is 
not  in  phase  over  the  horizontal  plane  but  has  a hexagonal  pattern 
and  produces  wave  trains  180°  out  of  phase  with  each  other  at  the 
center  of  the  hexagonal  cell  and  its  edge.  Both  wave  trains  develop 
into  shocks  and  the  horizontal  average  has  shocks  coming  twice  as 
fast  as  in  a single  wave  train.  The  fluctuating  component  of  the 
fluid  variables  develops  a square  wave  profile.  This  too  is  reason- 
able, because  when  the  square  wave  fluctuating  component  is  added 
with  appropriate  sign  to  the  saw  tooth  mean  component,  a saw 
tooth  shock  train  with  the  full  wave  length  results,  but  those 
shocks  at  the  cell  edge  are  displaced  a half  wave  length  in  height 
relative  to  those  at  the  cell  center.  On  the  other  hand,  we  have 
not  yet  been  able  to  verify  that  the  code  conserves  energy  properly. 
Further  work  on  finding  and  correcting  the  bugs  in  this  program 
is  being  carried  out. 

3.  Three  - Dimensional  Fluid  Dynamics 

Neither  one-dimensional  nor  even  "modal"  calculatio  .s  can 
Include  effects  of  horizontal  inhomogenelty  in  the  mean  structure 
of  the  solar  atmosphere.  To  do  that  at  least  two-and  preferably 
three-dimensional  calculations  are  needed.  Drs.  Auer,  Klein  and 
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myself  Degan  development  of  a three-dimensional  fluid  dynamic 
computer  code. 

Our  approach  is  to  represent  the  three-dimensional  difference 
operator  for  the  fluid  dynamic  equations  as  the  product  of  one- 
dimensional  operators.  Consider  for  example  the  continuity  equation 


a p a 

— + z . — ( P u ) = 0 

at  ax. 

1 

Let  represent  the  operator 

n+1/6  n At  a 

P = P ( p 

6 a x^ 


(10) 
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The  operators  for  different  directions  don't  commute  so  the  operator 
for  a full  time  step  is  taken  to  be 


L L L L L L (12) 

X y z z y X 

This  method  of  "time-splitting"  reduces  the  complicated  three- 
dimensional  differencing  problem  to  a product  of  simple  one-dimensional 
difference  algorithms.  An  additional  bonus  is  the  somewhat  larger  time 
step  permitted  by  the  Courant  condition. 

The  scheme  actually  used  to  transport  the  fluid  variables  in 
time  is  that  devised  by  MacCormack  (1971).  This  is  similar  to  the 
two-step  Lax-Wendroff  method  of  Rlchtmyer  and  Martin  (1969)  but 
(i)  always  uses  the  same  grid  points  (so  there  is  no  staggering 
in  space)  and  (ii)  changes  the  direction  of  spatial  differencing 
between  predictor  and  corrector,  so  one  is  always  upstream  and 
the  other  downstream.  This  Increases  the  dissipation  and  makes  the 
scheme  very  stable.  Additional  non-linear  stabilization  is  provided 
where  needed  by  the  Flux  Corrected  Transport  (FCT)  algorithm  of 
Boris  (Book,  Boris  and  Hain,  1975).  In  their  method  extra  diffusion 
is  added  to  the  transport  scheme  to  make  all  variables  vary 
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smoothly.  This  extra  diffusion  is  then  cancelled  out  everywhere 
except  where  that  would  lead  to  ripples.  The  process  of  limiting 
the  cancellation  where  ripples  would  arise  is  nonlinear.  We  are 
currently  developing  graphics  to  study  the  output  from  the  cal- 
culations and  a scheme  to  calculate  the  self -gravitational  attraction 
of  the  fluid. 

This  3D  code  will  be  tested  on  three  problems. 

(J)  Spherical  blast  wave  in  uniform  and  decreasing  density 
medium.  Comparsion  will  be  made  with  the  Sedov  similiarlty 
solutions . 

(ii)  Accretion  in  a rotating  disk  about  a point  mass  without 
self-gravity.  This  will  test  the  viscosity  and  angular  momen- 

I 

turn  transport  in  the  code.  Comparison  will  be  made  to  the  | 

i 

analytic  solutions  of  Lyndon-Bell  and  Pringle  (1974).  j 

! 

(iii)  Spherical  accretion  problem.  This  will  test  the  self 
gravity  scheme.  Comparison  will  be  made  to  the  similarity 
solutions  of  Larson  (1969). 
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B.  Waves  In  the  Transition  Region 


Skylab  observations  have  shovm  evidence  of  acoustic  pulses  In  the 
transition  region  (final  report  on  AF  contract  F19628-75-C-0013) , but 
there  Is  no  evidence  of  the  "f ive-mlnute"  oscillation  (Vernazza  et.  al. 
1975).  Also  there  Is  great  doubt  about  whether  acoustic  waves  supply 
sufficient  energy  to  heat  the  upper  chromosphere,  transition  region  and 
corona  (Jordan,  1977;  Athay  and  White,  1978;  Brunner  and  McWhlrter,  1978). 
Dr.  J.  V.«rnazza  and  I continued  our  study  of  vertically  pr^-’f; ^gating  acoustic 
waves  In  the  transition  region  In  order  to  Investigate  these  problems.  Our 
original  calculations  assumed  the  waves  were  Isothermal.  More  realistic 
calculations  have  now  been  performed. 

We  have  used  a one-dlmenslonal  Lagranglan  fluid  dynamic  code,  the 
energy  equation  Includes  electron  conduction  and  optically  thin  radiative 
losses  (Cox  and  Tucker,  1969).  The  Initial  structure  of  the  transition 
region  is  determined  primarily  by  the  condition  of  constant  conductive  flux. 
For  most  part  the  transition  region  is  very  finely  zoned  to  achieve 

4 

high  resolution.  The  one  exception  Is  the  region  between  4x10  °K  and 
3xl0^°K  which  Is  covered  by  only  one  zone  of  thickness  150  km.  This  will 
produce  some  unreal  numerical  wave  reflection,  but  the  effect  seems  to 
be  small.  The  waves  were  driven  by  a piston  1.5Mm  below  the  top  of  the 
convection  zone. 

Results  for  300s  waves  show: 

(I)  The  velocity  amplitude  is  nearly  constant  but  decreases  slightly 
from  the  middle  to  the  top  of  the  transition  region  (figure  1). 

(II)  The  density  amplitude  decreases  by  a factor  of  3 between  the 
top  of  the  chromosphere  and  the  transition  region  (figure  2). 
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Preliminary  results  on  acoustic  wave  fluxes  and  velocity  amplitudes 

confirm  Jordan's  results  that  acoustic  waves  can’t  supply  sufficient 

energy  to  heat  the  corona.  We  have  looked  at  wave  trains  with  75s, 

150s  and  300s  periods.  Waves,  with  velocity  amplitudes  less  than 

±1  km/s  at  the  temperature  minimum  have  energy  fluxes  less  than 
4 2 

2x10  erg/cm  s in  the  transition  region.  The  energy  flux  is 

A 9 

larger  for  150s  than  either  75s  or  300s  waves.  This  flux  is  2x10  erg/cm  s 

5 2 

too  small  to  supply  the  coronal  energy  needs  of  3x10  erg/cm  s. 

Another  interesting  result  is  that  the  velocity  amplitude  in  the 
transition  region  is  only  about  ± 20  km/s  when  it  is  ±1  km/s 
at  the  temperature  minimum.  This  transition  region  value  is  lower 
than  inferred  from  line  widths  by  Boland  et  al.  (1975)  even  though 
the  temperature  minimum  value  is  large  (Canfield  and  Beckers  1977). 
Furthermore,  the  velocity  amplitude  in  the  transition  region  is  very 
insensitive  to  the  driving  amplitude.  It  changed  by  less  than  25% 
for  a 300%  change  in  the  amplitude  at  the  temperature  minimum  and  the 
driving  piston. 
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C.  Minicomputer 
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A Nova  3/D  minicomputer  was  purchased  primarily  with  funds  from 
NSF,  but  with  a diskette  unit  purchased  with  funds  from  this  contract. 

It  was  delivered  toward  the  end  of  February  1978. 

When  operational  it  will  greatly  facilitate  work  under  this  contract: 

(a)  Time  dependent  calculations  will  be  run  for  very  long  times  - until 
statistical  steady  state  is  reached.  Such  calculations  are  prohibitively 
expensive  when  one  must  pay  by  the  hour  for  computer  time  (b;  A high 
level  editor  will  make  program  development  easier,  (c)  A phone  link 
to  the  CDC  7600  and  Cray  I computer  at  the  National  Center  for  Atmos- 
pheric Research  in  Boulder,  CO  will  be  implemented.  This  link  will 
enable  calculations  to  be  performed  with  these  super  high  speed  large 
memory  computers. 
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16. 


Figure  1:  Velocity  vs.  time  at  three  different  temperatures:  Solid  line  ^ 

4 i 

Is  for  10  "K  (below  the  transition  region) . Dashed  line  Is 

for  3.5  X 10^  “K  (middle  of  the  transition  region).  Dotted  j 

line  Is  for  9 x 10^  "K  (above  the  transition  region) . Waves  j 

are  driven  by  a piston  1500  km  below  the  photosphere  with  a 

-3  i 

velocity  amplitude  10  the  local  sound  speed.  ^ 

Figure  2:  Density  vs.  time  for  the  same  cases  as  In  figure  1. 
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